home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / Helix-turn-helix ƒ / gmd.txt < prev    next >
Encoding:
Text File  |  1994-04-01  |  25.8 KB  |  833 lines  |  [TEXT/R*ch]
  1.  
  2.                             DISTRIBUTION INFORMATION
  3.                               ON GENOMIC MAP DESIGN
  4.                                     (GMD)
  5.  
  6. Programs are now available to assist in the design of genome
  7. mapping experiments ("contig mapping"). The theory behind the GMD
  8. programs for designing "contig mapping" experiments can be found in
  9.  
  10. Fu, Y.-X., W.E. Timberlake, and J. Arnold (1992). "On the design
  11. of genome mapping experiments using short synthetic
  12. oligonucleotides", Biometrics, in press.
  13.  
  14. Any published use of these programs should cite this reference.
  15. This journal has a one year backlog, and it may be as late as 1993
  16. before the article appears.
  17.  
  18. OBTAINING THE SOFTWARE: The software is only distributed via
  19. Internet using EMAIL. Please send an EMAIL request to:
  20.  
  21.                     ARNOLD%GANDAL.DNET@SERVER.UGA.EDU
  22.  
  23. if you wish copies of the programs. I will EMAIL you three FORTRAN
  24. programs written to generate the tables in Fu et al. (1992). Each program
  25. is accompanied by a documentation file explaining the program's
  26. use.
  27.  
  28. USING THE SOFTWARE WITHOUT THE PROGRAMS: The programs also have been
  29. incorporated into a DNA sequence analysis package (Arnold et al., 1986),
  30. and can be accessed directly on the Biological Sequence/Structure
  31. Computational Facility (BS/SCF). Contact Dr. Weise for a guest account 
  32. at:
  33.                     WEISE%GANDAL.DNET@SERVER.UGA.EDU
  34.  
  35. OBTAINING FURTHER DOCUMENTATION: The best source of documentation
  36. is the paper by Fu et al. (1992). A reprint can be obtained by
  37. writing:
  38.                     Dr. Jonathan Arnold
  39.                     Genetics Department
  40.                     University of Georgia
  41.                     Athens, GA 30602
  42. or by emailing:
  43.                     ARNOLD%GANDAL.DNET@SERVER.UGA.EDU
  44.  
  45. SOFTWARE SUPPORT IN THE USE OF THE PROGRAMS: If you have questions about
  46. the programs, please contact Dr. Yun-xin Fu, who wrote the programs:
  47.  
  48.                     FU@GSBS18.GS.UTH.TMC.EDU
  49.  
  50. --------------------------------------------------------------------------------
  51.  
  52.  
  53.         DOCUMENTATION FOR TABLE1.EXE in DNASEQ
  54.  
  55.  
  56. This program provides information for designing genome mapping
  57. experiments. The map to be generated can be a RFLP map or
  58. physical map. The program generates Table 1 or variations
  59. on it in the reference:
  60.  
  61. Fu,  Y.X., W.E. Timberlake, and J. Arnold (1992). On the
  62. design of genome mapping experiments using short synthetic
  63. oligonucleotides, Biometrics, in press
  64.  
  65. This program is used to select library size or the number of
  66. RFLPs to map a chromosome. The program provides coverage statistics associated
  67. with each choice of library size (or number of RFLPs). The example
  68. is in terms of a physical mapping experiment.
  69.  
  70. To convert the example below into terms for an RFLP mapping experiment,
  71. make a proportional correspondence between the size of chromosome
  72. in map units and the target density of RFLPs (i.e. one per centimorgan)
  73. on the one hand with the size of the chromosome in kilobases and
  74. the insert size of a cloning vector. A probe is some short sequence
  75. motif that can be found at substantial frequency among clones in
  76. the library or associated with the collection of clones used to
  77. screen for the RFLPs. This short sequence motif may be detected
  78. by synthetic oligonucleotides or restriction enzymes, for example.
  79. For RFLP mapping take the chromosome size in map units
  80. and convert into kilobases (i.e. 2000 kb) and the insert
  81. size as (target density of RFLPs/chromosome size in map units)X2000.
  82.  
  83.  
  84. PROGRAM INPUT:
  85.  
  86. The program first asks you for the chromosome
  87. size and clone length (i.e. insert size or target density
  88. on RFLP map). The answer should be in kilobases, but if
  89. you are creating an RFLP map, please rescale proportionally
  90. the target density and length in map units so that the
  91. length in map units corresponds to kilobases.
  92.  
  93. Example: 400 40
  94.  
  95. The program then asks you for the clone range. That is, it is
  96. requesting the smallest library size, the largest library
  97. size, and increment in library size you would find
  98. acceptable. This question allows
  99. you to consider a range of values for the library size
  100. (number of clones) and to step through this range in
  101. steps of specified size.
  102.  
  103. Example: 40 50 5
  104.  
  105. PROGRAM OUTPUT:
  106.  
  107. The program then displays the predictions about the
  108. library on the terminal screen and also writes the results to a
  109. file with a name like FOR010.dat. The output from the
  110. example above is shown below and annotated:
  111.  
  112. Example:
  113.  
  114. The program repeats the input values for each library
  115. imagined. In the example, there were three libraries considered during
  116. input. The notation (cir) refers to a circular (bacterial)
  117. chromosome. The notation (Lin) refers to a linear
  118. chromosome. The notation (3a) refers to a circular
  119. chromosome; The notation (3b) refers to a linear
  120. chromosome. The notation (3c) refers to the interior
  121. of a linear chromosome. The notation (n.c) also
  122. refers to a circular chromosome, and the notation
  123. (c.l), to a linear chromosome.
  124.  
  125.  Genome size and clone length=       400        40
  126.             number of clones =        40
  127.            A: (cir) and (Lin)=   0.77858   0.77858     <----(expected
  128. proportion of clones overlapping the 3' end of a random clone)
  129.            B: (3a), (3b),(3c)=   0.985   0.988   0.944 <----(expected 
  130. proportion of coverage)
  131.            C: (n.c) and (n.l)=     3.9     4.1         <----(expected
  132. number of clones downstream of a random clone)
  133.            D: (c.c) and (c.l)=     1.0     1.5         <----(expected
  134. number of contigs)
  135.  
  136.  Genome size and clone length=       400        40
  137.             number of clones =        45
  138.            A: (cir) and (Lin)=   0.79962   0.79962
  139.            B: (3a), (3b),(3c)=   0.991   0.993   0.952
  140.            C: (n.c) and (n.l)=     4.4     4.6
  141.            D: (c.c) and (c.l)=     1.0     1.3
  142.  
  143.  Genome size and clone length=       400        40
  144.             number of clones =        50
  145.            A: (cir) and (Lin)=   0.81789   0.81789
  146.            B: (3a), (3b),(3c)=   0.995   0.996   0.958
  147.            C: (n.c) and (n.l)=     4.9     5.1
  148.            D: (c.c) and (c.l)=     1.0     1.2
  149.  
  150.  
  151.  
  152. EXITING PROGRAM:
  153.  
  154. Enter a string of zeros:
  155.  
  156. 0 0 0
  157.  
  158. Enter a second string of zeros:
  159.  
  160. 0 0 0
  161.  
  162. You should exit the program and return to the menu.
  163.  
  164.                     Jonathan Arnold
  165.  
  166. --------------------------------------------------------------------------------
  167.  
  168. c    Program used to generate Table 1
  169. c    Adapted from Nclone.for
  170. c
  171. c    Calculate the prob. of overlapping, En,
  172. c    E( Cn ) etc.
  173. c
  174.          program clone
  175.          implicit real*8 (a-h,o-z)
  176.          real*8  lark(-1:2000),fac(0:30000),lark2(-1:2000)
  177.          Real*8  g0,g1
  178.  
  179.          fac(0) = 0.0d0
  180.          Do k=1,30000,1
  181.              fac(k) = fac(k-1) + dlog( dble(k))
  182.          End Do         
  183.     Write(5,'(/A62/)') 'The results shown on screen are also
  184.      & in the file :for010.dat'
  185.          
  186.  10        print*,' Input genome size  and clone length (in kb)'
  187.          read(5,*) n,m
  188.          if( n .eq. 0 ) stop 
  189.  
  190.  20       print*,'Input clone range: start(say 10), end(say 50)
  191.      & and step (say 5)'
  192.          read(5,*) mk1,mk2,mk3
  193.          if( mk1 .eq. 0) go to 10
  194.          
  195.          do mkk = mk1,mk2,mk3
  196.  
  197.               mk = mkk-1
  198.              Ec0 = 1.0d0 - g0(N,M,mkk,2.0d0)
  199.              x1   =  g1(N,M,mkk+1,1.0d0)
  200.              x2   =  g1(N,M,mkk+1,2.0d0)
  201.              Ec1 = 1.0d0 - 2.0d0*(x1-x2)/(mkk+1)-x2
  202.              Ec2 = 1.0d0 - 2.0d0*(1.0d0-x2)/(mkk+1)-x2
  203.  
  204.              x0  = g0(N,M,1,2.0d0)
  205.  
  206.              xx0 = g1(N,M,2,2.0d0)
  207.              xx0 = (1.0d0-xx0)/2 + xx0
  208.  
  209.               xp = dlog(x0)
  210.              xq = dlog(1.0d0 - x0)
  211.  
  212.               xxp = dlog(xx0)
  213.              xxq = dlog(1.0d0 - xx0)
  214.  
  215.              en = 0.0d0
  216.              en1 = 0.0d0
  217.              Do k=0,mk,1
  218.                  x = fac(mk)-fac(k)-fac(mk-k)
  219.                  en = en + k*dexp( x + (mk-k)*xp + k*xq )
  220.                  en1 = en1 + k*dexp( x + (mk-k)*xxp + k*xxq )
  221.              end Do
  222.              Econtig0 = mkk*dexp( mk*xp )
  223.              Econtig1 = (mkk-1)*dexp( mk*xxp )+1
  224.  
  225.              lark(0) = g0(N,M,mkk-1,0.0d0)
  226.              lark2(0) = g1(N,M,mkk,0.0d0)
  227.              Do k=1,M,1
  228.                        x1 = 2.0d0*k/M
  229.                  xx = g0(N,M,mkk-1,x1)
  230.                  lark(k) =  xx
  231.  
  232.                  xx = g1(N,M,mkk,x1)
  233.                  lark2(k) = (1.0d0 -xx)/mkk + xx
  234.              end do
  235.  
  236.              yy =0.0d0
  237.              yy2=0.0d0
  238.  
  239.              ww =0.0d0
  240.              ww2=0.0d0
  241.              Do k=0,M-1,1
  242.                  ab = lark(k)-lark(k+1)
  243.                  ac = lark2(k)-lark2(k+1)
  244.                  ab2 = ab/(1.0d0-lark(M))
  245.                  ac2 = ab/(1.0d0-lark(M))
  246.  
  247.                  yy  =  yy + ab*k
  248.             yy2 = yy2 + ab2*k
  249.                  ww  =  ww + ac*k
  250.             ww2 = ww2 + ac2*k
  251.  
  252. c                 write(1,'(i6,4f10.5,i8)') k,ab,ab2,ac,ac2,mkk
  253. c                 write(5,'(i6,4f10.5,i8)') k,ab,ab2,ac,ac2,mkk
  254.  
  255.              end do
  256.  
  257.              v1= dmax1(Econtig0,1.0d0)
  258.              v2= dmax1(Econtig1,1.0d0)
  259.              Do lk=5,10,5
  260.              write(lk,'(/A30,2i10)') 'Genome size and clone length=',n,m
  261.              write(lk,'(A30,2i10)') 'number of clones =', mkk
  262.              write(lk,'(A30,2f10.5)') ' A: (cir) and (Lin)=',
  263.      +                    1.0d0-yy2/M,1-ww2/M
  264.              write(lk,'(a30,3f8.3)') '  B: (3a), (3b),(3c)=',Ec0,Ec1,Ec2
  265.              write(lk,'(a30,2f8.1)') '  C: (n.c) and (n.l)=',en,en1
  266.              write(lk,'(a30,2f8.1)') '  D: (c.c) and (c.l)=',v1,v2
  267.              end do
  268.  
  269.          end do
  270.  
  271.          go to 20
  272.          end
  273. c-------------------------------------------------------
  274. c    function 
  275. c------------------------------------------------------
  276.          Real*8 function g0(N,M,nn,x)
  277.          real*8 x,xx
  278.          integer nn,N,M
  279.  
  280.          xx = 1.0d0 - x*M/(2*N)
  281.          g0 = dexp( nn*dlog(xx) )
  282.  
  283.          end
  284. c------------------------------------------------------
  285.          Real*8 function g1(N,M,nn,x)
  286.          real*8 x,xx
  287.          integer nn,N,M
  288.  
  289.          xx = 1.0d0 - x*M/(2*(N-M))
  290.          g1 = dexp( nn*dlog(xx) )
  291.  
  292.          end
  293.  
  294. --------------------------------------------------------------------------------
  295.  
  296.  
  297.         DOCUMENTATION FOR TABLE2.EXE in DNASEQ
  298.  
  299.  
  300. This program provides information for designing genome mapping
  301. experiments. The map to be generated can be a RFLP map or
  302. physical map. The program generates Table 2 or variations
  303. on it in the reference:
  304.  
  305. Fu,  Y.X., W.E. Timberlake, and J. Arnold (1992). On the
  306. design of genome mapping experiments using short synthetic
  307. oligonucleotides, Biometrics, in press
  308.  
  309. This program is used to select the number of probes to
  310. be used to detect overlaps between clones. The program provides 
  311. the power to detect overlap as a function of the number
  312. of probes, the probability of clone/probe hybridization,
  313. the significance level (frequency of false positives tolerable).
  314.  
  315. To convert the example below into terms for an RFLP experiment,
  316. make a proportional correspondence between chromosome size in
  317. map units and the target density of RFLPs (i.e. one per centimorgan)
  318. on the one hand with chromosome size in kilobases and the insert
  319. size of a cloning vector on the other hand. A probe is some short sequence
  320. motif that can be found at substantial frequency among clones in
  321. the library or associated with clones used to screen for the RFLPs.
  322. This short sequence motif may be detected by synthetic oligonucleotides
  323. or restriction enzymes, for example. For RFLP mapping, take the
  324. chromosome size in map units and convert into kilobases
  325. (.i.e. 2000 kb) and the insert size, as 
  326. (target density of RFLPs in map units/chromosome size in map units)X2000.
  327.  
  328. PROGRAM INPUT:
  329.  
  330. The program first asks you for the insert size of your
  331. cloning vector (or target spacing on your RFLP map) and
  332. the probability of hybridization between clone and probe.
  333.  
  334. Example: 40 .40
  335.  
  336. The program then asks you for a range of probe numbers. That is, it is
  337. requesting the smallest number of probes to be used, the largest number
  338. of probes possible, and an increment in the number of probes. This question
  339. allows you to consider a range of values for the number of probes used
  340. and to step through this range in steps of specified size.
  341.  
  342. Example: 30 40 5
  343.  
  344. PROGRAM OUTPUT:
  345.  
  346. The program then writes the predictions about the
  347. experiment to a file with a name like FOR010.dat. The output from the
  348. example above is shown below and annotated:
  349.  
  350. Example:
  351.         Table 2, Powers of detecting two overlapping clones 
  352.                         Overlapping  (%)>=       0      20      50
  353.                         Probe length (bp)=       9
  354.  
  355.                           significance levels
  356.  
  357.         1.0d-7   1.0d-6   1.0d-5   1.0d-4   1.0d-3   0.0100   0.0250   0.0500
  358.  
  359.                         Clone length (kb)=      40
  360.                         prob. of hybrid. =   0.4000
  361.    30   0.1096   0.1096   0.1668   0.2865   0.3491   0.4803   0.5487   0.6181
  362.         0.1371   0.1371   0.2085   0.3580   0.4359   0.5960   0.6753   0.7506
  363.         0.2192   0.2192   0.3330   0.5626   0.6710   0.8476   0.9079   0.9490
  364.    35   0.1425   0.1924   0.2436   0.2962   0.4057   0.5213   0.5811   0.6416
  365.         0.1782   0.2405   0.3045   0.3702   0.5063   0.6459   0.7141   0.7784
  366.         0.2849   0.3841   0.4842   0.5834   0.7642   0.8955   0.9376   0.9656
  367.    40   0.1677   0.2120   0.2572   0.3510   0.4495   0.5527   0.6059   0.6596
  368.         0.2096   0.2649   0.3215   0.4387   0.5606   0.6841   0.7437   0.7997
  369.         0.3352   0.4230   0.5116   0.6835   0.8290   0.9269   0.9568   0.9762
  370.  
  371. There are a triple of numbers associated with each combination. Consider
  372. the last column, where the significance level (the acceptable level
  373. to you of false positives). The number .6181 is the power (chance) of detecting
  374. any overlap when the number of probes used is 30. The number
  375. .7506 is the power (chance) of detecting an overlap of at
  376. least 20%. The last number .9490 is the power (chance) of detecting
  377. an overlap of at least 50%. These three detection probabilities, power,
  378. are for an experiment with a probability of clone/probe hybridization
  379. of .4 and insert size of 40kb.
  380.  
  381. EXITING PROGRAM:
  382.  
  383. Enter a string of zeros:
  384.  
  385. 0 0 0
  386.  
  387. Enter a second string of zeros:
  388.  
  389. 0 0 0
  390.  
  391. You should exit the program and return to the menu.
  392.  
  393.                     Jonathan Arnold
  394.  
  395. --------------------------------------------------------------------------------
  396.  
  397. c    adapted from semi.for to compute table 2.
  398. c    difference is that here three overlapping 0,20,50 are
  399. c    combined to give one table.
  400. c    march 27,90
  401. c
  402.          program semi
  403.          implicit real*8 (a-h),(o-z)
  404.          integer crit(0:8),crit2(0:8)
  405.          real*8 fac(0:5000),qq(0:3000),siglevel(0:8),power(0:8),
  406.      &  power2(0:8)
  407.  
  408.          write(5,'(/a55/)') 'Output is stored in the file :For011.dat'
  409.  
  410.          fac(0) = 0.0d0
  411.          do k=1,5000,1
  412.              fac(k) = dlog( dble(k)) + fac(k-1)
  413.          end do         
  414.          lpr = 9
  415.  
  416.          do k=11,11,4
  417.          write(k,'(/a60)') 'Table 2, Powers of detecting two overlapping
  418.      & clones '
  419.          write(k,'(a42,3i8)')    ' Overlapping  (%)>=', 0,20,50
  420.          write(k,'(a42,i8)')    ' Probe length (bp)=', lpr
  421.          write(k,'(/a45/)')    ' significance levels'
  422.          write(k,'(5x,8a9/)') '1.0d-7','1.0d-6','1.0d-5','1.0d-4' ,
  423.      &'1.0d-3',
  424.      +                '0.0100','0.0250','0.0500'
  425.          end do
  426.  
  427. 10         Print*,' Input clone length(in kb) and hybri prob.'
  428.          Read(*,*) lcl,ph
  429.          if( lcl .eq. 0) stop
  430.          Print*,' Input the range of probe numbers: start, end and step'
  431.          read(*,*) mk1,mk2,mk3
  432.  
  433. c         do lll = 30,40,5
  434.              lll=lcl
  435.              lcl=lll*1000
  436.              w = lcl- lpr+1
  437.               write(11,'(a42,i8)')    ' Clone length (kb)=', lcl/1000
  438. c             do ph=0.30d0,0.51d0,0.05d0
  439.               write(11,'(a42,f9.4)')  ' prob. of hybrid. =', ph
  440.  
  441.              pq = 2.0d0*ph*(1.0d0-ph)
  442.              pp = 1.0d0 -pq
  443.              a = dlog( 1.0d0-ph )/w
  444.  
  445.          Do nprobe = mk1,mk2,mk3
  446.               do lmn = 1,3,1
  447.          
  448.              if(lmn .eq.1 ) mper = 0
  449.              if(lmn .eq.2 ) mper = 20*(lcl/100)
  450.              if(lmn .eq.3 ) mper = 50*(lcl/100)
  451.          
  452.              kp = lcl-125
  453.              num = 0
  454.              d2  = dlog(2.0d0)
  455.              Do k=0,nprobe,1
  456.                  qq(k)=0.0d0
  457.               end do
  458.  
  459.              Do loc = 1,500,1
  460.                  if( kp .lt. mper ) go to 100
  461.                  k  = kp-lpr+1
  462.                   q = d2+ a*w + dlog( 1.0d0 - dexp(a*(w-k)) )
  463.                  p = dlog(1.0d0 - dexp(q) )
  464.                  Do i=0,nprobe,1
  465.                       qq(i) = qq(i)+ 
  466.      +                    dexp( dble(i)*q + dble(nprobe-i)*p )
  467.                  end do
  468.                  kp = kp-250
  469.              end do
  470.  
  471.  100    num = loc-1
  472.  
  473. c         print*,' the number of location considered =',num
  474.          siglevel(0) = 0.050d0
  475.          siglevel(1) = 0.025d0
  476.          siglevel(2) = 0.010d0
  477.        Do k=0,7,1
  478.              power2(k) = 0.0d0
  479.              power(k) = 0.0d0
  480.              crit(k)  = 0
  481.              crit2(k)  = 0
  482.              if( k.gt.2) siglevel(k)=siglevel(k-1)*0.1d0
  483.          end do
  484.          cy = 0.0d0
  485.          cx = 0.0d0
  486.          p = dlog(pp)
  487.          q = dlog(pq)
  488.  
  489.     Do k=0,nprobe,1
  490.              ww =  fac(nprobe)-fac(k)-fac(nprobe-k)
  491.              qq(k) = qq(k)/num
  492.              xx = dexp(ww)*qq(k)
  493.              yy = dexp( ww + k*q + (nprobe-k)*p )
  494.              cx = cx + xx
  495.              cy = cy + yy
  496.  
  497.              Do i=0,7,1
  498.                  if( cy .lt. siglevel(i) ) then
  499.                      power(i) = cx
  500.                      crit(i)  = k
  501.                   end if
  502.                  if((1.0d0-cx .lt. siglevel(i)).and.(crit2(i).eq.0))then
  503.                      power2(i) = 1.0d0-cy
  504.                      crit2(i)  = k
  505.                   end if
  506.              end do
  507.  
  508. c             Do i=10,10,5
  509. c                    write(unit=i,fmt='(i4,2f10.6,i5,f8.4,2i9)') k,
  510. c     +            xx,yy,nprobe,ph,lcl,lpr
  511. c             end do
  512.          end do
  513. c         write(5,'(i5,8f9.4)') nprobe,(power(k),k=7,0,-1)
  514. c         write(5,'(5x,8f9.4)')        (power2(k),k=7,0,-1)
  515.          if( lmn .eq.1) then
  516.               write(11,'(i5,8f9.4)') nprobe,(power(k),k=7,0,-1)
  517.          else     
  518.         write(11,'(5x,8f9.4)')        (power(k),k=7,0,-1)
  519.          end if
  520. c         write(2,'(i5,8i9)') nprobe, (crit(k),k=7,0,-1)
  521.  
  522.         end do
  523.          end do
  524.  
  525. c         end do
  526. c         end do
  527.          go to 10
  528.  
  529.          end
  530.  
  531. --------------------------------------------------------------------------------
  532.  
  533.  
  534.         DOCUMENTATION FOR TABLE3.EXE in DNASEQ
  535.  
  536.  
  537. This program provides information for designing genome mapping
  538. experiments. The map to be generated can be a RFLP map or
  539. physical map. The program generates Table 3 or variations
  540. on it in the reference:
  541.  
  542. Fu,  Y.X., W.E. Timberlake, and J. Arnold (1992). On the
  543. design of genome mapping experiments using short synthetic
  544. oligonucleotides, Biometrics, in press
  545.  
  546. This program is used to select the number of probes and clones to
  547. be used to detect overlaps between clones in order to
  548. enlarge a contig. The program provides 
  549. the power to detect overlap between a given contig
  550. and an overlapping clone downstream for varying
  551. significance levels, the number
  552. of probes, number of clones, and chromosome size.
  553.  
  554. To convert the example below into terms for an RFLP experiment
  555. make a proportional correspondence between chromosome size
  556. in map units and the target desnity of RFLPs (i.e. one per
  557. centimorgan). on the one hand with chromosome size in kilobases
  558. and the insert size of a cloning vector on the other hand.
  559. A probe is some short sequence motif that can be found at
  560. substantial frequency among clones in the library or associated
  561. with the collection of clones used to screen for the RFLPs.
  562. This short sequence motif may be detected by restriction enzymes
  563. or synthetic oligonucleotides, for example. For RFLP mapping,
  564. take the chromosome size  in map units and convert
  565. into kilobases (i.e 2000 kb) and the insert size as
  566. (target density of RFLPs in map units/chromosome size in map units)X2000.
  567.  
  568. PROGRAM INPUT:
  569.  
  570. The program first asks you for the chromosome size, clone length (insert
  571. size of cloning vector), and number of clones in the library.
  572.  
  573. Example: 400 40 50
  574.  
  575. The next question is a request for the probability of
  576. clone/probe hybridization:
  577.  
  578. Example: .40
  579.  
  580. The final question is for a range of probe numbers. That is, it is
  581. requesting the smallest number of probes to be used, the largest number
  582. of probes possible, and an increment in the number of probes. This question
  583. allows you to consider a range of values for the number of probes used
  584. and to step through this range in steps of specified size.
  585.  
  586. Example: 30 40 5
  587.  
  588. PROGRAM OUTPUT:
  589.  
  590.   The program then writes the predictions about the
  591. experiment to a file with a name like FOR012.dat. The output from the
  592. example above is shown below and annotated:
  593.  
  594. Example:
  595.  
  596.      Table 3. powers of detecting the nearest right-hand side overlapping clone
  597.                     With  Probe length   =       9  bp
  598.  
  599.                           significance levels
  600.  
  601.         1.0d-7   1.0d-6   1.0d-5   1.0d-4   1.0d-3   0.0100   0.0250   0.0500
  602.  
  603.                         prob.of hybri.   =   0.400
  604.                         Genome size      =     400  kb
  605.                         Clone length     =      40  kb
  606.                         Clone number     =      50
  607.    30   0.3830   0.3830   0.5212   0.7195   0.7887   0.8846   0.9167   0.9410
  608.    35   0.4718   0.5784   0.6659   0.7373   0.8417   0.9087   0.9320   0.9503
  609.    40   0.5328   0.6186   0.6905   0.8000   0.8746   0.9243   0.9423   0.9566
  610.  
  611. There is now a single number associated with each combination. Consider
  612. the last row. The number .9410 is the power (chance) of detecting
  613. any overlap downstream of a contig when the number of probes used is 30
  614. and the significance level (frequency of false positives) is .05. 
  615. The rest of the characteristics of the experiment are listed 
  616. in the Table legend, i.e. the probability of hybridization,
  617. chromosome size, clone length, and clone number (library size).
  618.  
  619. EXITING PROGRAM:
  620.  
  621. Enter a string of zeros:
  622.  
  623. 0 0 0
  624.  
  625. Enter a second string of zeros:
  626.  
  627. 0 0 0
  628.  
  629. You should exit the program and return to the menu.
  630.                         
  631.                     Jonathan Arnold
  632.  
  633. --------------------------------------------------------------------------------
  634.  
  635. C***************************************************************************
  636. c    adapted from npower.for to produce table3 in the genome 
  637. c    mapping paper.
  638. C    March 27,1990.
  639. C    Changed for Janothan on June 20.
  640. c***************************************************************************
  641.          program power
  642.          implicit real*8 (a-h),(o-z)
  643.          integer crit(0:8),crit2(0:8)
  644.          real*8 fac(0:1000),qq(0:300),siglevel(0:8),pow(0:8),
  645.      &  pow2(0:8)
  646.          real*8  lark(-1:2000),lark2(-1:2000),density(0:1000)
  647.          real*8  g0,g1
  648.  
  649.          write(*,'(/a30/)') 'The output is in for012.dat'
  650.  
  651.          kcode = 0
  652. c    kcode = 1 will be linear genome
  653.  
  654.         fac(0) = 0.0d0
  655.          do k=1,1000,1
  656.              fac(k) = dlog( dble(k)) + fac(k-1)
  657.          end do         
  658.          lpr = 9
  659. c         ph  = 0.40d0
  660.   
  661.          write(12,'(/a79)') 'Table 3. powers of detecting the nearest
  662.      & right-hand side overlapping clone'
  663.          write(12,'(a42,i8,a4)')    ' With  Probe length   =', lpr,' bp'
  664.          write(12,'(/a45/)')    ' significance levels'
  665.          write(12,'(5x,8a9/)') '1.0d-7','1.0d-6','1.0d-5','1.0d-4' ,
  666.      &'1.0d-3',
  667.      +                '0.0100','0.0250','0.0500'
  668. c-------------------------------------------------------------------------
  669. 10         Print*,' Input genome size (kb),clone length(kb) and
  670.      & clone number'
  671.          read(*,*) n,m,mm
  672.  
  673.          if( n .eq. 0 ) stop
  674.            print*,' input prob.of hybri.'
  675.          read(*,*) ph
  676.  
  677.            print*,' input probe number range: start, end and step'
  678.          read(*,*) mk1,mk2,mk3
  679.  
  680. c         Do nleng = 1,3,1
  681. c             if(nleng .eq. 1) n = 300
  682. c             if(nleng .eq. 2) n = 2000
  683. c             if(nleng .eq. 3) n = 10000
  684.               write(12,'(a42,f8.3)')    ' prob.of hybri.   =', ph
  685.               write(12,'(a42,i8,a4)')    ' Genome size      =', n,' kb'
  686.                   n = n*2
  687. c
  688. c             do mleng=1,3,1
  689. c             if( mleng.eq.1) m = 30
  690. c             if( mleng.eq.2) m = 35
  691. c             if( mleng.eq.3) m = 40
  692.  
  693.               write(12,'(a42,i8,a4)')    ' Clone length     =', m,' kb'
  694.               lcl= m*1000
  695.              m = m*2
  696.          
  697. c             do mmleng=1,4,1
  698. c             if( mmleng .eq.1) mm = 50
  699. c             if( mmleng .eq.2) mm = 100
  700. c             if( mmleng .eq.3) mm = 200
  701. c             if( mmleng .eq.4) mm = 400
  702.               write(12,'(a42,i8)')    ' Clone number     =', mm
  703.  
  704. c
  705. c    to obtain the distance distribution
  706. c
  707.          w = lcl- lpr+1
  708.          pq = 2.0d0*ph*(1.0d0-ph)
  709.          pp = 1.0d0 -pq
  710.          a = dlog( 1.0d0-ph )/w
  711.  
  712.  
  713.          if( kcode .eq. 0) then
  714.               lark(0) = g0(N,M,mm-1,0.0d0)
  715.              yy  = g0(N,M,mm-1,2.0d0)
  716.          else     
  717.         lark(0) = g1(N,M,mm,0.0d0)
  718.              yy  = g1(N,M,mm,2.0d0)
  719.              yy  = (1.0d0-yy)/(mm+1.0d0)+ yy 
  720.          end if
  721.          Do k=1,M,1
  722.                    x1 = 2.0d0*k/M
  723.              if( kcode .eq.0) then
  724.                   xx = g0(N,M,mm-1,x1)
  725.                  lark(k) =  xx
  726.         else        
  727.                       xx = g1(N,M,mm,x1)
  728.                   lark(k) = (1.0d0 -xx)/mm + xx
  729.              end if     
  730.              density(k-1) = (lark(k-1)-lark(k))/(1.0d0-yy)
  731.     end do
  732. c
  733. c    computing powers
  734. c
  735.  
  736.          Do nprobe = mk1,mk2,mk3
  737.              kp = lcl-250
  738.              num = 0
  739.              d2  = dlog(2.0d0)
  740.              Do k=0,nprobe,1
  741.                  qq(k)=0.0d0
  742.               end do
  743.  
  744.              x00 = 0.0d0
  745.              Do loc = 1,500,1
  746.                  if( kp .lt. 0 ) go to 100
  747.                  x00 = x00 + density(loc-1)
  748. c                 print*, loc,kp,density(loc-1),x00
  749.  
  750.                  k  = kp-lpr+1
  751.                   q = d2+ a*w + dlog( 1.0d0 - dexp(a*(w-k)) )
  752.                  p = dlog(1.0d0 - dexp(q) )
  753.                  Do i=0,nprobe,1
  754.                       qq(i) = qq(i)+ density(loc-1)* 
  755.      +                    dexp( dble(i)*q + dble(nprobe-i)*p )
  756.                  end do
  757.                  kp = kp-500
  758.              end do
  759.  
  760.  100    num = loc-1
  761.  
  762. c         print*,' the number of location considered =',num
  763.          siglevel(0) = 0.050d0
  764.          siglevel(1) = 0.025d0
  765.          siglevel(2) = 0.010d0
  766.        Do k=0,7,1
  767.              pow2(k) = 0.0d0
  768.              pow(k) = 0.0d0
  769.              crit(k)  = 0
  770.              crit2(k)  = 0
  771.              if( k.gt.2) siglevel(k)=siglevel(k-1)*0.1d0
  772.          end do
  773.          cy = 0.0d0
  774.          cx = 0.0d0
  775.          p = dlog(pp)
  776.          q = dlog(pq)
  777.  
  778.     Do k=0,nprobe,1
  779.              ww =  fac(nprobe)-fac(k)-fac(nprobe-k)
  780.              xx = dexp(ww)*qq(k)
  781.              yy = dexp( ww + k*q + (nprobe-k)*p )
  782.              cx = cx + xx
  783.              cy = cy + yy
  784.  
  785.              Do i=0,7,1
  786.                  if( cy .lt. siglevel(i) ) then
  787.                      pow(i) = cx
  788.                      crit(i)  = k
  789.                   end if
  790.                  if((1.0d0-cx .lt. siglevel(i)).and.(crit2(i).eq.0))then
  791.                      pow2(i) = 1.0d0-cy
  792.                      crit2(i)  = k
  793.                   end if
  794.              end do
  795.  
  796. c             Do i=10,10,5
  797. c                    write(unit=i,fmt='(i4,2f10.6,i5,f8.4,4i9)') k,
  798. c     +            xx,yy,nprobe,ph,n,m,mm,lpr
  799. c             end do
  800.          end do
  801. c         write(5,'(i5,8f9.4)') nprobe,(pow(k),k=7,0,-1)
  802. c         write(5,'(5x,8f9.4)')        (pow2(k),k=7,0,-1)
  803.          write(12,'(i5,8f9.4)') nprobe,(pow(k),k=7,0,-1)
  804. c         write(12,'(5x,8f9.4)')        (pow2(k),k=7,0,-1)
  805. c         write(2,'(i5,8i9)')   nprobe,(crit(k),k=7,0,-1)
  806.          end do
  807.  
  808. c         end do
  809. c         end do
  810. c         end do
  811.          go to 10
  812.          end
  813. c-------------------------------------------------------
  814. c    function 
  815. c------------------------------------------------------
  816.          Real*8 function g0(N,M,nn,x)
  817.          real*8 x,xx
  818.          integer nn,N,M
  819.  
  820.          xx = 1.0d0 - x*M/(2*N)
  821.          g0 = dexp( nn*dlog(xx) )
  822.  
  823.          end
  824. c------------------------------------------------------
  825.          Real*8 function g1(N,M,nn,x)
  826.          real*8 x,xx
  827.          integer nn,N,M
  828.  
  829.          xx = 1.0d0 - x*M/(2*(N-M))
  830.          g1 = dexp( nn*dlog(xx) )
  831.  
  832.          end
  833.